NetBoxR Tutorial

Overview

The netboxr package composes a number of functions to retrive and process genetic data from large-scale genomics projects (e.g. TCGA projects) including from mutations, copy number alterations, gene expression and DNA methylation. The netboxr package implements NetBox algorithm in R package. NetBox algorithm integrates genetic alterations with literature-curated pathway knowledge to identify pathway modules in cancer. NetBox algorithm uses (1) global network null model and (2) local network null model to access the statistic significance of the discovered pathway modules.

Basics

Installation

library(devtools)
install_bitbucket(repo="mil2041/netboxr", ref="master", build_vignette=TRUE)

Getting Started

Load netboxr package:

library(netboxr)

A list of all accessible vignettes and methods is available with the following command:

help.search("netboxr")

For help on any netboxr package functions, use one of the following command formats:

help(geneConnector)
?geneConnector

Example of Cerami et al. PLoS Comp Bio 2010

This is an example to reproduce the network discovered on Cerami et al.(2010).

The results presented here are comparable to the those from Cerami et al. 2010 though the unadjusted p-values for linker genes are not the same. It is because the unadjusted p-value of linker genes in Cerami et al. 2010 were calculated by the probabiliy of the observed data point, Pr(X). The netboxr used the probability of an observed or more extreme assuming the null hypothesis is true, Pr(X>=x|H), as unadjusted p-value for linker genes. The final number of linker genes after FDR correction are the same between netboxr result and original Cerami et al. 2010.

Load Human Interactions Network (HIN) network

Load pre-defined HIN network and simplify the interactions by removing loops and duplicated interactions in the network. The netowork after reduction contains 9264 nodes and 68111 interactions.

data(netbox2010)
sifNetwork <- netbox2010$network
graphReduced <- networkSimplify(sifNetwork,directed = FALSE)      

Loading network of 9264 nodes and 157780 interactions
Treated as undirected network 
Removing multiple interactions and loops
Returning network of 9264 nodes and 68111 interactions

Load altered gene list

The altered gene list contains 517 candidates from mutations and copy number alterations.

geneList <- as.character(netbox2010$geneList) 
length(geneList)

[1] 517

Map altered gene list on HIN network

The geneConnector function in the netboxr package takes altered gene list as input and maps the genes on the curated network to find the local processes represented by the gene list.

## Use Benjamini-Hochberg method to do multiple hypothesis 
## correction for linker candidates.
## Use edge-betweeness method to detect community structure in the network. 
threshold <- 0.05
results <- geneConnector(geneList=geneList,
                        networkGraph=graphReduced,
                        directed=FALSE,
                        pValueAdj="BH",
                        pValueCutoff=threshold,
                        communityMethod="ebc",
                        keepIsolatedNodes=FALSE)

274 / 517 candidate nodes match the name in the network of 9264 nodes 
Only test neighbor nodes with local degree equals or exceeds 2
Multiple hypothesis corrections for 892 neighbor nodes in the network
For p-value 0.05 cut-off, 6 nodes were included as linker nodes
Connecting 274 candidate nodes and 6 linker nodes
Remove 208 isolated candidate nodes from the input
Final network contains 72 nodes and 152 interactions
Detecting modules using "edge betweeness" method

# Check the p-value of the selected linker
linkerDF <- results$neighborData
linkerDF[linkerDF$pValueFDR<threshold,]

##
## The geneConnector function returns a list of data frames. 
names(results)

[1] "netboxGraph"      "netboxCommunity"  "netboxOutput"     "nodeType"        
[5] "moduleMembership" "neighborData"    

# plot graph with the Fruchterman-Reingold layout algorithm
plot(results$netboxCommunity,results$netboxGraph, layout=layout_with_fr) 

Consistency with Previously Published Results

The GBM result by netboxr identified exactly the same linker genes (6 linker genes), the same number of modules (10 modules) and the same genes in each identified module as GBM result in Cerami et al. 2010.

The results of netboxr are consistent with previous implementation of the NetBox algorithm. The RB1 and PIK3R1 modules are clearly represented in the figure. For example, the RB1 module contains genes in blue color and enclosed by light orange circle. The PIK3R1 module contains genes in orange color and enclosed by pink circle.

Statistical Significance of Discovered Modules

NetBox algorithm used (1) global network null model and (2) local network null model to access the statistic significance of the discovered modules.

Global Network Null Model

The global network null model calculates the empirical p-value as the number of times (over a set of iterations) the size of the largest module in the network coming from the same number of randomly selected genes (number of genes is 274 in this example) equals or exceeds the size of the largest module in the observed network. Suggested iterations is 1000.

## This function will need a lot of time to complete. 
globalTest <- globalNullModel(netboxGraph=results$netboxGraph, networkGraph=graphReduced, iterations=10, numOfGenes = 274)

Local Network Null Model

Local network null model evaluates the deviation of modularity in the observed network from modularity distribution in the random network. For each interaction, a random network is produced from local re-wiring of literature curated network. It means all nodes in the network kept the same degree of connections but connect to new neighbors randomly. Suggested iterations is 1000.

localTest <- localNullModel(netboxGraph=results$netboxGraph, iterations=1000)

###########
Based on 1000 random trails
Random networks: mean modularity = 0.295
Random networks: sd modularity = 0.059
Observed network modularity is: 0.519
Observed network modularity z-score is: 3.804
One-tail p-value is: 7.131e-05

Through 1000 iterations, we can obtain the mean and the standard deviation of modularity in the local network null model. Using the mean (~0.3) and the standard deviation (0.06), we can covert the observed modularity in the network (0.519) into a Z-score (~3.8). From the Z-score, we can calculate one-tail p-value. If one-tail pvalue is less than 0.05, the observed modularity is significantly different from random. In the histogram, the blue region is the distribution of modularity in the local network null model. The red vertical line is the observed modularity in the NetBox results.

h<-hist(localTest$randomModularityScore,breaks=35,plot=FALSE)
h$density = h$counts/sum(h$counts)
plot(h,freq=FALSE,ylim=c(0,0.1),xlim=c(0.1,0.6), col="lightblue")
abline(v=localTest$modularityScoreObs,col="red")

The global null model is used to assess the global connectivity (number of nodes and edges) of the largest module in the identified network compared with the same number but randomly selected gene list. The local null model is used to assess the network modularity in the identified network compared with random re-wired network.

View Module Membership

The table below shows the module memberships for all genes.

DT::datatable(results$moduleMembership, rownames = FALSE)

Write NetBox Output to Files

## Write results for further visilaztion in the cytoscape software. 
##
## network.sif file is the NetBox algorithm output in SIF format.  
write.table(results$netboxOutput, file="network.sif", sep="\t", quote=FALSE, col.names=FALSE, row.names=FALSE)
##
## netighborList.txt file contains the information of all neighbor nodes. 
write.table(results$neighborData, file="neighborList.txt", sep="\t", quote=FALSE, col.names=TRUE, row.names=FALSE)
##
## community.membership.txt file indicates the identified pathway module numbers.
write.table(results$moduleMembership, file="community.membership.txt", sep="\t", quote=FALSE, col.names=FALSE, row.names=FALSE)
##
## nodeType.txt file indicates the node is "linker" node or "candidate" node. 
write.table(results$nodeType,file="nodeType.txt", sep="\t", quote=FALSE, col.names=FALSE, row.names=FALSE)

Term Enrichment in Modules using Gene Ontology (GO) Analysis

After module identification, one main task is understanding the biological processes that may be represented by the returned modules. Here we use the Bioncoductor clusterProfiler to do an enrichment analysis using GO Biological Process terms on a selected module.

library(clusterProfiler)

Loading required package: DOSE

DOSE v3.4.0  For help: https://guangchuangyu.github.io/DOSE

If you use DOSE in published research, please cite:
Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609

clusterProfiler v3.6.0  For help: https://guangchuangyu.github.io/clusterProfiler

If you use clusterProfiler in published research, please cite:
Guangchuang Yu., Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.

Attaching package: ‘clusterProfiler’

The following object is masked from ‘package:igraph’:

    simplify

library(org.Hs.eg.db)
module <- 6
selectedModule <- results$moduleMembership[results$moduleMembership$membership == module,]
geneList <-selectedModule$geneSymbol
# Check available ID types in for the org.Hs.eg.db annotation package
keytypes(org.Hs.eg.db)

 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"
 [6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"    
[11] "GO"           "GOALL"        "IPI"          "MAP"          "OMIM"        
[16] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
[21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"       "UNIGENE"     
[26] "UNIPROT"     

ids <- bitr(geneList, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb="org.Hs.eg.db")

'select()' returned 1:1 mapping between keys and columns
25% of input gene IDs are fail to map...

head(ids)

ego <- enrichGO(gene = ids$ENTREZID,
                OrgDb = org.Hs.eg.db,
                ont = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
                readable = TRUE)

Enrichment Results

head(ego)

Visualize Enrichment Results

dotplot(ego)

Alternative Module Discovery Methods

In netboxr, we used the Girvan-Newman algorithm (communityMethod=“ebc”) as the default method to detect community membership in the identified network. The Girvan-Newman algorithm iteratativly removes the edge in the network with highest edge betweeness until no edges left. When the identified network contains many edges, the Girvan-Newman algorithm will spend a large amount of time to remove edges and re-calucalte the edge betweenese score in the network. If the user cannot get the community detection result in reasonable time, we suggest to switch to leading eigenvector method (communityMethod=“lec”) for community detection. Users can check original papers of the Girvan-Newman algorithm and leading eigenvector method for more details.

Alternative Pathway Data using PaxtoolsR

Users can load alternative pathway data from the Pathway Commons repository using the paxtoolsr package from Bioconductor. This pathway data represents an update to the Pathway Commons data used in the original 2010 NetBox publication. Below is an example that makes use of data from the Reactome pathway database.

NOTE: Downloaded data is automatically cached to avoid unnecessary downloads.

library(paxtoolsr)
filename <- "PathwayCommons.8.reactome.EXTENDED_BINARY_SIF.hgnc.txt.gz"
sif <- downloadPc2(filename, version="8")
# NOTE: Run without filename to see a list of available files
#sif <- downloadPc2()
# Filter interactions for specific types
interactionTypes <- getSifInteractionCategories()
filteredSif <- filterSif(sif$edges, interactionTypes=interactionTypes[["BetweenProteins"]])
filteredSif <- filteredSif[(filteredSif$INTERACTION_TYPE %in% "in-complex-with"), ]
# Re-run NetBox algorithm with new network
graphReduced <- networkSimplify(filteredSif, directed=FALSE)      

Loading network of 6440 nodes and 105767 interactions
Treated as undirected network 
Removing multiple interactions and loops
Returning network of 6440 nodes and 105767 interactions

geneList <- as.character(netbox2010$geneList) 
threshold <- 0.05
pcResults <- geneConnector(geneList=geneList,
                           networkGraph=graphReduced,
                           directed=FALSE,
                           pValueAdj="BH",
                           pValueCutoff=threshold,
                           communityMethod="lec",
                           keepIsolatedNodes=FALSE)

181 / 517 candidate nodes match the name in the network of 6440 nodes 
Only test neighbor nodes with local degree equals or exceeds 2
Multiple hypothesis corrections for 1132 neighbor nodes in the network
For p-value 0.05 cut-off, 31 nodes were included as linker nodes
Connecting 181 candidate nodes and 31 linker nodes
Remove 90 isolated candidate nodes from the input
Final network contains 122 nodes and 423 interactions
Detecting modules using "leading eigenvector" method

# Check the p-value of the selected linker
linkerDF <- results$neighborData
linkerDF[linkerDF$pValueFDR<threshold,]

# The geneConnector function returns a list of data frames. 
names(results)

[1] "netboxGraph"      "netboxCommunity"  "netboxOutput"     "nodeType"        
[5] "moduleMembership" "neighborData"    

# plot graph with the Fruchterman-Reingold layout algorithm
plot(results$netboxCommunity,results$netboxGraph, layout=layout_with_fr) 

References

• Cerami E, Demir E, Schultz N, Taylor BS, Sander C (2010) Automated Network Analysis Identifies Core Pathways in Glioblastoma. PLoS ONE 5(2): e8918. doi:10.1371/journal.pone.0008918
• Cerami EG, Gross BE, Demir E, Rodchenkov I, Babur O, Anwar N, Schultz N, Bader GD, Sander C. Pathway Commons, a web resource for biological pathway data. Nucleic Acids Res. 2011 Jan;39(Database issue):D685-90. doi:10.1093/nar/gkq1039. Epub 2010 Nov 10.

Session Information

sessionInfo()

R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
 [1] parallel  stats4    grid      stats     graphics  grDevices utils     datasets 
 [9] methods   base     

other attached packages:
 [1] clusterProfiler_3.6.0 DOSE_3.4.0            netboxr_0.99.4       
 [4] org.Hs.eg.db_3.5.0    AnnotationDbi_1.40.0  IRanges_2.12.0       
 [7] S4Vectors_0.16.0      Biobase_2.38.0        BiocGenerics_0.24.0  
[10] ComplexHeatmap_1.17.1 paxtoolsr_1.6.7       XML_3.98-1.9         
[13] rJava_0.9-9           igraph_1.1.2          R.utils_2.6.0        
[16] R.oo_1.21.0           R.methodsS3_1.7.1     BiocStyle_2.6.0      
[19] knitr_1.20           

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.16         circlize_0.4.2       tidyr_0.7.2          GO.db_3.5.0         
 [5] rprojroot_1.2        digest_0.6.12        R6_2.2.2             plyr_1.8.4          
 [9] backports_1.1.1      RSQLite_2.0          evaluate_0.10.1      httr_1.3.1          
[13] ggplot2_2.2.1        GlobalOptions_0.0.12 rlang_0.1.4          lazyeval_0.2.1      
[17] data.table_1.10.4-3  blob_1.1.0           GetoptLong_0.1.6     DT_0.2              
[21] qvalue_2.10.0        rmarkdown_1.8        labeling_0.3         splines_3.4.2       
[25] BiocParallel_1.12.0  stringr_1.3.0        htmlwidgets_0.9      bit_1.1-12          
[29] munsell_0.4.3        fgsea_1.4.0          compiler_3.4.2       pkgconfig_2.0.1     
[33] base64enc_0.1-3      shape_1.4.3          htmltools_0.3.6      tibble_1.3.4        
[37] gridExtra_2.3        jsonlite_1.5         gtable_0.2.0         DBI_0.7             
[41] magrittr_1.5         scales_0.5.0         stringi_1.1.7        GOSemSim_2.4.0      
[45] reshape2_1.4.2       DO.db_2.9            rvcheck_0.0.9        fastmatch_1.1-0     
[49] rjson_0.2.15         RColorBrewer_1.1-2   tools_3.4.2          bit64_0.9-7         
[53] glue_1.2.0           purrr_0.2.4          yaml_2.1.18          colorspace_1.3-2    
[57] memoise_1.1.0